
function fext2(fname,ext)
character*80 fname,fext2
character*3 ext
parameter(n=80)

fext2=fname

do 10 i=30,1,-1
  if(fext2(i:i).eq.'.') then
    fext2(i:n)=char(0)
    go to 11
 endif
10 continue
11 continue

do 20 i=30,1,-1
  if(fext2(i:i).ne.' ' .and. fext2(i:i).ne.char(0)) then
    fext2(i+1:i+4)='.'//ext
    go to 21
  endif
20 continue
   write(*,*) 'erro em fext'

21 continue
   return
end
                 
subroutine readgrd_lake(nunit,var,i0,i1,j0,j1)
!
! reads 2d array from a grd file
!
implicit real*8(a-h,o-z)
dimension var(i0:i1,j0:j1)
character*4 dsaa

read(nunit,'(a4)') dsaa
read(nunit,*) nx,ny
read(nunit,*) xmin,xmax
read(nunit,*) ymin,ymax
read(nunit,*) zmin,zmax
      
do j=j0,j1
  read(nunit,*) (var(i,j),i=i0,i1)
enddo
      
return
end
      
subroutine wrigrd_lake(nunit,z,dx,dy,i0,i1,j0,j1)

implicit none
integer i0,i1,j0,j1,nunit,i,j
real(8) z(i0:i1,j0:j1)
real(8) xmin,xmax,ymin,ymax,zmin,zmax,dx,dy
      
zmin=z(i0,j0)
zmax=z(i0,j0)

do j=j0,j1
  do i=i0,i1
    zmin=min(zmin,z(i,j))
    zmax=max(zmax,z(i,j))
  enddo
enddo

xmin=float(i0)*dx
xmax=float(i1)*dx
ymin=float(j0)*dy
ymax=float(j1)*dy
      
write(nunit,'(a4)') 'DSAA'
write(nunit,*) (i1-i0+1),(j1-j0+1)
write(nunit,*) xmin,xmax
write(nunit,*) ymin,ymax
write(nunit,*) zmin,zmax

do j=j0,j1
  write(nunit,*) (z(i,j),i=i0,i1)
enddo
     
return
end


SUBROUTINE MON_OUT(ix,iy,nx,ny,year,month,day,hour)
use ARRAYS
use INOUT_PARAMETERS, only : &
& lake_mon_out_unit_min, &
& lake_mon_out_unit_max
implicit none

! Input variables
! Reals
real(8), intent(in) :: hour

! Integers
integer(4), intent(in) :: year
integer(4), intent(in) :: month
integer(4), intent(in) :: day
integer(4), intent(in) :: ix
integer(4), intent(in) :: iy
integer(4), intent(in) :: nx
integer(4), intent(in) :: ny

! Local variables
! Reals
real(8), allocatable :: accum_var (:,:,:,:)
real(8), allocatable :: var       (:,:,:,:)
real(8), allocatable :: tsteps    (:,:)

real(8), external:: DZETA

! Integers
integer(4), allocatable :: month_old(:,:)
integer(4) :: out_unit = lake_mon_out_unit_min
integer(4) :: i ! Loop index

! Characters
character :: month1*2
character :: year1*4
character :: day1*2
character :: hour1*2
character :: outpath*60
character :: timestring*6
character :: coords_point*6

! Logicals
logical :: firstcall

common /out/ outpath
data firstcall /.true./

SAVE

if (firstcall) then
  allocate (month_old(1:nx, 1:ny))
  allocate (tsteps    (1:nx, 1:ny))  
  allocate (var      (1:4, 1:nx, 1:ny, 1:M+1) )
  allocate (accum_var(1:4, 1:nx, 1:ny, 1:M+1) )
  month_old(:,:) = month
  tsteps   (:,:) = 0.d0  
  accum_var(:,:,:,:) = 0.d0
  var      (:,:,:,:) = 0.d0
endif

var(1,ix,iy,1:M+1) = Tw1 (1:M+1)
var(2,ix,iy,1:M+1) = Sal1(1:M+1)
var(3,ix,iy,1:M)   = E1  (1:M)
var(4,ix,iy,1:M)   = eps1(1:M)

accum_var(:,ix,iy,:) = accum_var(:,ix,iy,:) + var(:,ix,iy,:)
tsteps(ix,iy) = tsteps(ix,iy) + 1.d0

if (month_old(ix,iy)/=month) then
  accum_var(:,ix,iy,:) = accum_var(:,ix,iy,:)/tsteps(ix,iy)
  call DATEMINUS (1,year,month,day,hour,year1,month1,day1,hour1)
  call TIMESTR(6,year1,month1,day1,hour1,timestring)
  call CHECK_UNIT(lake_mon_out_unit_min,lake_mon_out_unit_max,out_unit)
  write (coords_point,'(2i3)') ix, iy
  open(out_unit,file=outpath(1:len_trim(outpath))//'monthly/'// &
 & 'Profiles'//coords_point//timestring//'.dat', status='unknown')
  write (out_unit,*) '1 - depth, m' 
  write (out_unit,*) '2 - temperature, C'
  write (out_unit,*) '3 - salinity, kg/kg' 
  write (out_unit,*) '4 - turbulent kinetic energy, m**2/s**2'
  write (out_unit,*) '5 - disspation rate, m**2/s**3'
  do i=1,M+1 
    write (out_unit,'(f10.6,f9.3,3e12.4)') &
    & -dzeta_int(i)*h1, accum_var(1:4,ix,iy,i)
  enddo
  close(out_unit)
  month_old(ix,iy) = month
  tsteps(ix,iy) = 0.d0
  accum_var(:,ix,iy,:) = 0.d0
endif

if (firstcall) firstcall=.false.
END SUBROUTINE MON_OUT


SUBROUTINE DAY_OUT(ix,iy,nx,ny,year,month,day,hour)
use ARRAYS
use INOUT_PARAMETERS, only : &
& lake_day_out_unit_min, &
& lake_day_out_unit_max
implicit none

! Input variables
! Reals
real(8), intent(in) :: hour

! Integers
integer(4), intent(in) :: year
integer(4), intent(in) :: month
integer(4), intent(in) :: day
integer(4), intent(in) :: ix
integer(4), intent(in) :: iy
integer(4), intent(in) :: nx
integer(4), intent(in) :: ny

! Local variables
! Reals
real(8), allocatable :: tsteps   (:,:)
real(8), allocatable :: accum_var(:,:,:,:)
real(8), allocatable :: var      (:,:,:,:)

real(8), external :: DZETA

! Integers
integer(4), allocatable :: day_old(:,:)
integer(4) :: i ! Loop index
integer(4) :: out_unit = lake_day_out_unit_min


character :: month1*2
character :: year1*4
character :: day1*2
character :: hour1*2
character :: outpath*60
character :: timestring*8
character :: coords_point*6

logical :: firstcall

common /out/ outpath
data firstcall /.true./
      
SAVE

if (firstcall) then
  allocate (day_old   (1:nx, 1:ny))
  allocate (tsteps    (1:nx, 1:ny))  
  allocate (var       (1:4, 1:nx, 1:ny, 1:M+1) )
  allocate (accum_var (1:4, 1:nx, 1:ny, 1:M+1) )
  day_old  (:,:) = day
  tsteps   (:,:) = 0.d0  
  accum_var(:,:,:,:) = 0.d0
  var      (:,:,:,:) = 0.d0
endif

var(1,ix,iy,1:M+1) = Tw1 (1:M+1)
var(2,ix,iy,1:M+1) = Sal1(1:M+1)
var(3,ix,iy,1:M)   = E1  (1:M)
var(4,ix,iy,1:M)   = eps1(1:M)

accum_var(:,ix,iy,:) = accum_var(:,ix,iy,:) + var(:,ix,iy,:)
tsteps(ix,iy) = tsteps(ix,iy) + 1.d0

if (day_old(ix,iy)/=day) then
  accum_var(:,ix,iy,:) = accum_var(:,ix,iy,:)/tsteps(ix,iy)
  call DATEMINUS (2,year,month,day,hour,year1,month1,day1,hour1)
  call TIMESTR (8,year1,month1,day1,hour1,timestring)
  call CHECK_UNIT(lake_day_out_unit_min,lake_day_out_unit_max,out_unit)
  write (coords_point, '(2i3)') ix, iy
  open(out_unit, file=outpath(1:len_trim(outpath))//'daily/'// &
  & 'Profiles'//coords_point//timestring//'.dat',status='unknown')
  write (out_unit,*) '1 - depth, m' 
  write (out_unit,*) '2 - temperature, C'
  write (out_unit,*) '3 - salinity, kg/kg' 
  write (out_unit,*) '4 - turbulent kinetic energy, m**2/s**2'
  write (out_unit,*) '5 - disspation rate, m**2/s**3'
  do i=1,M+1 
    write (out_unit,'(f10.6,f9.3,3e12.4)') &
    &  -dzeta_int(i)*h1, accum_var(1:4,ix,iy,i)
  enddo
  close(out_unit)
  day_old(ix,iy) = day
  accum_var(:,ix,iy,:) = 0.d0
  tsteps(ix,iy) = 0.d0
endif

if (firstcall) firstcall=.false.
END SUBROUTINE DAY_OUT


SUBROUTINE HOUR_OUT(ix,iy,nx,ny,year,month,day,hour)

use ARRAYS, only:  &
 & Tw1, Sal1,      &
 & E1 , eps1,      &
 & h1 ,            &
 & PEMF,           &
 & PT_DOWN, PSAL_DOWN, &
 & PDENS_DOWN,     &
 & row,            &
 & k_turb_T_flux,  &
 & T_massflux,     & 
 & H_mixed_layer,  & 
 & w_conv_scale,   &
 & T_conv_scale,   &
 & Buoyancy0,      &
 & S,              &
 & Gen,            &
 & TKE_turb_trans, &
 & k3_mid,         &
 & u1, v1,         &
 & S_integr_positive, &
 & S_integr_negative, &
 & Gen_integr, &
 & eps_integr, &
 & TKE_balance, &
 & H_entrainment, &
 & dzeta_int, dzeta_05int

use DRIVING_PARAMS, only : &
& M, &
& scale_output

use INOUT_PARAMETERS, only : &
& lake_hour_out_unit_min, &
& lake_hour_out_unit_max

implicit none

! Input variables
! Reals
real(8), intent(in) :: hour

! Integers
integer(4), intent(in) :: ix
integer(4), intent(in) :: iy
integer(4), intent(in) :: nx
integer(4), intent(in) :: ny
integer(4), intent(in) :: year
integer(4), intent(in) :: month
integer(4), intent(in) :: day

! Local variables
! Reals
real(8), allocatable :: accum_var (:,:,:,:)
real(8), allocatable :: var       (:,:,:,:)
real(8), allocatable :: var_scalar(:,:,:)
real(8), allocatable :: accum_var_scalar(:,:,:)

real(8), external :: DZETA

! Integers
integer(4), parameter :: n_var = 21
integer(4), parameter :: n_var_scalar = 6
integer(4), allocatable :: hour_old(:,:)
integer(4), allocatable :: tsteps(:,:)
integer(4), allocatable :: n_unit(:,:)
integer(4) :: i
integer(4) :: out_unit = lake_hour_out_unit_min

character :: month1*2
character :: year1*4
character :: day1*2
character :: hour1*2
character :: outpath*60
character :: timestring*10
character :: coords_point*6
character :: format_char*100

logical :: firstcall = .true.
logical, allocatable :: firstcall_ixiy(:,:)

common /out/ outpath
      
SAVE

if (firstcall) then
  allocate (accum_var (n_var,1:nx,1:ny,1:M+1) )
  allocate (var       (n_var,1:nx,1:ny,1:M+1) )
  allocate (accum_var_scalar (n_var_scalar,1:nx,1:ny) )
  allocate (var_scalar       (n_var_scalar,1:nx,1:ny) )  
  allocate (tsteps    (1:nx,1:ny) )
  allocate (hour_old  (1:nx,1:ny) )
  allocate (n_unit(1:nx,1:ny))
  allocate (firstcall_ixiy(1:nx,1:ny))
  accum_var(:,:,:,:) = 0.d0
  accum_var_scalar(:,:,:) = 0.d0
  var(:,:,:,:) = 0.d0
  var_scalar(:,:,:) = 0.d0
  tsteps(:,:) = 0
  hour_old(:,:) = int(hour)
  firstcall_ixiy(:,:) = .true.
  n_unit(:,:) = out_unit + 1
endif

var(1,ix,iy,1:M)   = PEMF          (1:M)
var(2,ix,iy,1:M+1) = PDENS_DOWN    (1:M+1)
var(3,ix,iy,1:M+1) = PT_DOWN       (1:M+1)
var(4,ix,iy,1:M+1) = PSAL_DOWN     (1:M+1)
var(5,ix,iy,1:M)   = k_turb_T_flux (1:M)
var(6,ix,iy,1:M)   = T_massflux    (1:M)
var(7,ix,iy,1:M+1) = row           (1:M+1)
var(8,ix,iy,1:M+1) = Tw1           (1:M+1)
var(9,ix,iy,1:M+1) = Sal1          (1:M+1)
var(10,ix,iy,1:M)  = E1            (1:M)
var(11,ix,iy,1:M)  = eps1          (1:M)

if (scale_output == 1) then
  var(12,ix,iy,1)    = H_mixed_layer ! Scale
  var(13,ix,iy,1)    = Buoyancy0     ! Scale
  var(14,ix,iy,1)    = w_conv_scale  ! Scale
  var(15,ix,iy,1)    = T_conv_scale  ! Scale
elseif (scale_output == 0) then
  var(12,ix,iy,1)    = 1.d0 ! No scaling
  var(13,ix,iy,1)    = 1.d0 ! No scaling
  var(14,ix,iy,1)    = 1.d0 ! No scaling
  var(15,ix,iy,1)    = 1.d0 ! No scaling
else
  print*, 'Scale_output must be 0 or 1: STOP'  
  STOP
endif
      
var(16,ix,iy,1:M)  = S              (1:M)
var(17,ix,iy,1:M)  = Gen            (1:M)
var(18,ix,iy,1:M)  = TKE_turb_trans (1:M)
      
var(19,ix,iy,2:M)  = k3_mid         (2:M)

var(20,ix,iy,1:M+1)  = u1(1:M+1)
var(21,ix,iy,1:M+1)  = v1(1:M+1)

var_scalar(1,ix,iy) = S_integr_positive
var_scalar(2,ix,iy) = S_integr_negative
var_scalar(3,ix,iy) = Gen_integr
var_scalar(4,ix,iy) = eps_integr
var_scalar(5,ix,iy) = TKE_balance
var_scalar(6,ix,iy) = H_entrainment

accum_var_scalar(:,ix,iy) = accum_var_scalar(:,ix,iy) + var_scalar(:,ix,iy)
accum_var(:,ix,iy,:) = accum_var(:,ix,iy,:) + var(:,ix,iy,:)
tsteps(ix,iy) = tsteps(ix,iy) + 1

if (int(hour)/=hour_old(ix,iy)) then
  accum_var(:,ix,iy,:) = accum_var(:,ix,iy,:) / real(tsteps(ix,iy))
  accum_var_scalar(:,ix,iy) = accum_var_scalar(:,ix,iy) / real(tsteps(ix,iy))
  call DATEMINUS(3,year,month,day,hour,year1,month1,day1,hour1)
  call TIMESTR(10,year1,month1,day1,hour1,timestring)
  call CHECK_UNIT(lake_hour_out_unit_min,lake_hour_out_unit_max,out_unit)
  write (coords_point, '(2i3)') ix, iy
  
! Writing to the file Profiles<yyyymmddhh>.dat  
  open(out_unit, file=outpath(1:len_trim(outpath))// &
  & 'hourly/'//'Profiles'//coords_point//timestring//'.dat', &
  &  status='unknown')
  write (out_unit,*) '1 - depth, normalized' 
  write (out_unit,*) '2 - temperature, normalized'
  write (out_unit,*) '3 - salinity, kg/kg' 
  write (out_unit,*) '4 - water density, kg/m**3'  
  write (out_unit,*) '5 - turbulent kinetic energy, normalized'
  write (out_unit,*) '6 - disspation rate, normalized'
  write (out_unit,*) '7 - eddy diffusivity (TKE**2/dissipation), m**2/s'
  write (out_unit,*) '8  - mass flux, m/s'
  write (out_unit,*) '9  - downdraft temperature, C' 
  write (out_unit,*) '10 - downdraft salinity, kg/kg' 
  write (out_unit,*) '11 - downdraft density, kg/m**3'
  write (out_unit,*) '12 - x-component of speed, m/s'
  write (out_unit,*) '13 - y-component of speed, m/s'
  write (out_unit,'(11(i6,4x))') 1,2,3,4,5,6,7,8,9,10,11 
  do i=1,M 
    write (out_unit,'(f14.6,f12.5,e10.3,f10.3,3e10.3,4e14.7,2f11.5)')      &
    & -dzeta_int(i)*h1                         /accum_var(12,ix,iy,1),  &
    & (accum_var(8,ix,iy,i) - scale_output*maxval(accum_var(8,ix,iy,:)) ) &
    &  /accum_var(15,ix,iy,1), &
    & accum_var(9,ix,iy,i),                                               &
    & accum_var(7,ix,iy,i),                                               &
    & accum_var(10,ix,iy,i)    / (accum_var(13,ix,iy,1)*accum_var(12,ix,iy,1) /       &
    &  accum_var(14,ix,iy,1) ),                                           &
    & accum_var(11,ix,iy,i)    / accum_var(13,ix,iy,1),                   &
    & accum_var(19,ix,iy,i),                                              & !accum_var(10,i)**2 / accum_var(11,i),
    & accum_var(1,ix,iy,i),                                               &
    & accum_var(3,ix,iy,i),                                               &
    & accum_var(4,ix,iy,i),                                               &
    & accum_var(2,ix,iy,i),                                               &
    & accum_var(20,ix,iy,i),                                              &
    & accum_var(21,ix,iy,i)
  enddo 
  close (out_unit)
  
! Writing to the file EDMF_profiles<yyyymmddhh>.dat
  open(out_unit, file=outpath(1:len_trim(outpath))//  &
  & 'hourly/'//'EDMF_profiles'//coords_point//timestring//'.dat',   &
  & status='unknown')
  write (out_unit,*) '1 - depth, normalized by mixed layer depth' 
  write (out_unit,*) '2 - "k - flux" of temperature, normalized'
  write (out_unit,*) '3 - mass flux of temperature,  normalized'
  write (out_unit,*) '4 - total flux of temperature, normalized'
  write (out_unit,'( 4(i10,6x) )') 1,2,3,4
  do i=1, M
    write (out_unit,'(f18.6,3e16.7)')                   &
    & -dzeta_05int(i)*h1/ accum_var(12,ix,iy,1),          &
    & accum_var(5,ix,iy,i)/(accum_var(14,ix,iy,1)*accum_var(15,ix,iy,1)), &
    & accum_var(6,ix,iy,i)/(accum_var(14,ix,iy,1)*accum_var(15,ix,iy,1)), &
    & (accum_var(5,ix,iy,i) + accum_var(6,ix,iy,i) ) /              &
    & (accum_var(14,ix,iy,1)*accum_var(15,ix,iy,1))
  enddo
  close (out_unit)
  
! Writing to the file TKE_budget<yyyymmddhh>.dat
  open(out_unit, file=outpath(1:len_trim(outpath))//    &
  & 'hourly/'//'TKE_budget'//coords_point//timestring//'.dat',        &
  &  status='unknown')
  write (out_unit,*) '1 - depth, normalised with mixed layer depth'
  write (out_unit,*) '2 - shear production, normalised with surface buoyancy flux'
  write (out_unit,*) '3 - buoyancy source,  normalised with surface buoyancy flux'
  write (out_unit,*) '4 - dissipation rate, normalised with surface buoyancy flux'
  write (out_unit,*) '5 - turbulent transport, normalised with surface buoyancy flux'
  write (out_unit,'( 5(i10,6x) )') 1,2,3,4,5
  do i = 1, M
    write (out_unit,'(f18.6,4e16.7)')          &
    & -dzeta_05int(i)*h1/ accum_var(12,ix,iy,1), &
    & accum_var(17,ix,iy,i) / accum_var(13,ix,iy,1), &
    & accum_var(16,ix,iy,i) / accum_var(13,ix,iy,1), &
    & accum_var(11,ix,iy,i) / accum_var(13,ix,iy,1), &
    & accum_var(18,ix,iy,i) / accum_var(13,ix,iy,1)
  enddo
!  write (out_unit, '(a)')
  close (out_unit)
  
  if (firstcall_ixiy(ix,iy)) then
    call CHECK_UNIT(lake_hour_out_unit_min,lake_hour_out_unit_max,n_unit(ix,iy))
    open(n_unit(ix,iy), file=outpath(1:len_trim(outpath))//           &
    & 'hourly/'//'TKE_integr_conv'//coords_point//timestring//'.dat', &
    &  status='unknown')
    write (n_unit(ix,iy),*) '1 - time,          hours'
    write (n_unit(ix,iy),*) '2 - H_entrainment, m    '
    write (n_unit(ix,iy),*) '3 - B_integr+,   m**3/s**3'
    write (n_unit(ix,iy),*) '4 - B_integr-,   m**3/s**3'
    write (n_unit(ix,iy),*) '5 - S_integr,    m**3/s**3'
    write (n_unit(ix,iy),*) '6 - eps_integr,  m**3/s**3'
    write (n_unit(ix,iy),*) '7 - TKE_balance, m**3/s**3'
    firstcall_ixiy(ix,iy) = .false.
  endif
  
  format_char = '(2f10.2, 5e15.6)'
  write (n_unit(ix,iy), format_char)  hour,       &
  &                    accum_var_scalar(6,ix,iy), & 
  &                    accum_var_scalar(1,ix,iy), &
  &                    accum_var_scalar(2,ix,iy), &
  &                    accum_var_scalar(3,ix,iy), &
  &                    accum_var_scalar(4,ix,iy), &
  &                    accum_var_scalar(5,ix,iy)
  
  hour_old(ix,iy) = int(hour)
  accum_var(:,ix,iy,:) = 0.d0
  accum_var_scalar(:,ix,iy) = 0.d0
  tsteps(ix,iy)  = 0
   
endif

if (firstcall) firstcall = .false.
!if (firstcall_ixiy(ix,iy)) firstcall_ixiy(ix,iy) = .false.
END SUBROUTINE HOUR_OUT          
      


SUBROUTINE EVERYSTEP_OUT(ix,iy,nx,ny)

use ARRAYS, only : &
& Tw1, &
& nstep, &
& M, &
& h1, &
& Sal1, &
& E1, &
& eps1, &
& k_turb_T_flux, &
& dzeta_int

use ARRAYS, only : &
 & S_integr_positive, &
 & S_integr_negative, &
 & Gen_integr, &
 & eps_integr, &
 & TKE_balance, &
 & E_integr, &
 & H_entrainment, &
 & time

use DRIVING_PARAMS, only : &
& everystep

use INOUT_PARAMETERS, only : &
& lake_everystep_out_unit_min, &
& lake_everystep_out_unit_max

implicit none

! Input variables
! Integers
integer(4), intent(in) :: ix
integer(4), intent(in) :: iy
integer(4), intent(in) :: nx
integer(4), intent(in) :: ny

! Local variables
! Reals
real(8), external :: DZETA
real(8), parameter :: ACC = 1.d-20
real(8), parameter :: hour_sec = 60.*60.

! Integers
integer(4), allocatable :: n_unit(:,:,:)
integer(4) :: i ! Loop index

! Characters
character :: outpath*60
character :: coords_point*6
character :: format_char*100

! Logicals
logical, allocatable :: firstcall(:,:)

common /out/ outpath
      
SAVE

if (.not.allocated(firstcall)) then
  allocate (firstcall(1:nx, 1:ny) )
  allocate (n_unit   (1:2, 1:nx, 1:ny) )
  firstcall = .true.
  n_unit = lake_everystep_out_unit_min
endif  
      
if (firstcall(ix,iy)) then
  write (coords_point, '(2i3)') ix, iy
  
  if (everystep /= 2) then
    call CHECK_UNIT(lake_everystep_out_unit_min,lake_everystep_out_unit_max, &
    & n_unit(1,ix,iy))
    open (n_unit(1,ix,iy), file=outpath(1:len_trim(outpath))//'everystep/'// &
    & 'Profiles'//coords_point//'.dat',  status='unknown')
    write (n_unit(1,ix,iy),*) '1 - depth, m' 
    write (n_unit(1,ix,iy),*) '2 - temperature, C'
    write (n_unit(1,ix,iy),*) '3 - salinity, kg/kg' 
    write (n_unit(1,ix,iy),*) '4 - turbulent kinetic energy, m**2/s**2'
    write (n_unit(1,ix,iy),*) '5 - disspation rate, m**2/s**3'
    write (n_unit(1,ix,iy),*) '6 - eddy diffusivity (TKE**2/dissipation), m**2/s'
    write (n_unit(1,ix,iy),*) '7 - k-flux of temperature, K*m/s'
  endif
    
  call CHECK_UNIT(lake_everystep_out_unit_min,lake_everystep_out_unit_max, &
  & n_unit(2,ix,iy))
  open (n_unit(2,ix,iy), file=outpath(1:len_trim(outpath))//'everystep/'// &
  & 'TKE_integr'//coords_point//'.dat',  status='unknown')
  write (n_unit(2,ix,iy),*) '1 - timestep              '
  write (n_unit(2,ix,iy),*) '2 - time,        hours    '
  write (n_unit(2,ix,iy),*) '3 - B_integr+,   m**3/s**3'
  write (n_unit(2,ix,iy),*) '4 - B_integr-,   m**3/s**3'
  write (n_unit(2,ix,iy),*) '5 - S_integr,    m**3/s**3'
  write (n_unit(2,ix,iy),*) '6 - eps_integr,  m**3/s**3'
  write (n_unit(2,ix,iy),*) '7 - TKE_balance, m**3/s**3'
  write (n_unit(2,ix,iy),*) '8 - E_integr,    m**3/s**2'
endif

if (everystep /= 2) then      
  write (n_unit(1,ix,iy),*) 'nstep = ', nstep
  format_char = '(f10.6,f9.3,5e12.4)'
  do i=1,M 
    write (n_unit(1,ix,iy), format_char) &
    & -dzeta_int(i)*h1,Tw1(i),Sal1(i),E1(i),eps1(i),  &
    & E1(i)**2/(eps1(i)+ACC), k_turb_T_flux(i)
  enddo
endif

format_char = '(i10, 2f10.2, 6e15.6)'
write(n_unit(2,ix,iy), format_char) &
& nstep, time/hour_sec, H_entrainment, S_integr_positive, S_integr_negative, &
& Gen_integr, eps_integr, TKE_balance, E_integr
      
if (firstcall(ix,iy)) firstcall(ix,iy)=.false.
END SUBROUTINE EVERYSTEP_OUT
      
      
SUBROUTINE SERIES_OUT(ix,iy,nx,ny,year,month,day,hour,tsw)

use DRIVING_PARAMS, only : &
& M, ns, &
& dt_out
use NUMERIC_PARAMS , only : &
& ms
use ARRAYS, only:  &
& row,             &
& Tw1,Tskin,       &
& Ti1,             &
& zinv,            &
& time,            &
& h1,l1,hs1,ls1,   &
& H_mixed_layer,   &
& H_entrainment,   &
& qwater, qsoil,   &
!& qwater2, qsoil2, & ! two-meth
& oxyg,            &
& febul, fplant,   &
& fdiff_lake_surf, &
& fdiff,           &
!& febul2,           & ! two-meth
!& fdiff_lake_surf2, & ! two-meth
!& fdiff2,           & ! two-meth
& rprod_total_newC, &
& rprod_total_oldC, &
& h_talik,          &
& T_0dim
use ATMOS,  only:    &
& eflux0_kinem,      &
& eflux0,            &
& turb_density_flux, &
& hw, xlew, tau,     &
& surfrad
use INOUT_PARAMETERS, only : &
& lake_series_out_unit_min, &
& lake_series_out_unit_max
use METH_OXYG_CONSTANTS, only : &
& molmass_ch4

implicit none

! Input variables
! Reals
real(8), intent(in) :: tsw
real(8), intent(in) :: hour

! Integers
integer(4), intent(in) :: ix, iy
integer(4), intent(in) :: nx, ny
integer(4), intent(in) :: year, month, day

! Local variables
! Reals
real(8), parameter :: hour_sec = 60.*60., day_sec = 24.*60.*60.
real(8), parameter :: mfs = molmass_ch4*8.64d+7 ! the transform multiplier for methane flux
                                                ! from mol/(m**2*s) to mg/(m**2*day)
real(8), parameter :: small_number = 1.d-5

real(8) :: T_mean
real(8) :: roughness,emissivity,albedo,aM,bM,relhums

real(8), external :: VARMEAN

! Integers
integer(4), parameter :: nfiles = 5
integer(4), allocatable :: n_unit(:,:,:)
integer(4), allocatable :: count_out(:,:)

! Characters
character :: outpath*60
character :: coords_point*6
character :: format_char*100

! Logicals
logical, allocatable :: firstcall(:,:)

common /out/ outpath
common /surface/ roughness,emissivity,albedo,aM,bM,relhums

real(8) :: dz(ms)
integer(4) :: itop
common /SOILDAT/ dz,itop
real(8) :: Tsn(1:ms), cs(1:ms)
common /SNOW_CHAR/ Tsn,cs
      
SAVE

if (.not.allocated(firstcall)) then
  allocate (firstcall(1:nx, 1:ny) )
  allocate (count_out(1:nx, 1:ny) )
  allocate (n_unit   (1:nfiles, 1:nx, 1:ny) )
  firstcall(:,:) = .true.
  n_unit = lake_series_out_unit_min
endif  
            
if (firstcall(ix,iy)) then
  write (coords_point, '(2i3)') ix, iy

  call CHECK_UNIT(lake_series_out_unit_min,lake_series_out_unit_max, &
  & n_unit(1,ix,iy))
  open (n_unit(1,ix,iy),file=outpath(1:len_trim(outpath))//'time_series/'// &
  & 'layers'//coords_point//'.dat',  status='unknown')
  write (n_unit(1,ix,iy),*)'Col. 1 - year'
  write (n_unit(1,ix,iy),*)'Col. 2 - month'
  write (n_unit(1,ix,iy),*)'Col. 3 - day'
  write (n_unit(1,ix,iy),*)'Col. 4 - hour'
  write (n_unit(1,ix,iy),*)'Col. 5 - the time from the start of integration, hours'
  write (n_unit(1,ix,iy),*)'Col. 6 - water layer thickness, m'
  write (n_unit(1,ix,iy),*)'Col. 7 - ice layer thickness,   m'
  write (n_unit(1,ix,iy),*)'Col. 8 - snow layer thickness,  m'
  write (n_unit(1,ix,iy),*)'Col. 9 - bottom ice thickness,  m'
    
  call CHECK_UNIT(lake_series_out_unit_min,lake_series_out_unit_max, &
  & n_unit(2,ix,iy))
  open (n_unit(2,ix,iy),file=outpath(1:len_trim(outpath))//'time_series/'// &
  & 'T_fluxes'//coords_point//'.dat',status='unknown')
  write (n_unit(2,ix,iy),*)'Col. 1 - year'
  write (n_unit(2,ix,iy),*)'Col. 2 - month'
  write (n_unit(2,ix,iy),*)'Col. 3 - day'
  write (n_unit(2,ix,iy),*)'Col. 4 - hour'
  write (n_unit(2,ix,iy),*)'Col. 5 - the time from the start of integration, hours'
  write (n_unit(2,ix,iy),*)'Col. 6 - surface temperature, C'
  write (n_unit(2,ix,iy),*)'Col. 7 - water skin temperature, C'
  write (n_unit(2,ix,iy),*)'Col. 8 - water surface temperature, C'
  write (n_unit(2,ix,iy),*)'Col. 9 - water bottom temperature, C'
  write (n_unit(2,ix,iy),*)'Col. 10 - mean temperature of water coloumn, C'
  write (n_unit(2,ix,iy),*)'Col. 11 - zero-dimensional model temperature, C'  
  write (n_unit(2,ix,iy),*)'Col. 12 - upper ice surface temperature, C'
  write (n_unit(2,ix,iy),*)'Col. 13 - upper snow surface temperature, C'
  write (n_unit(2,ix,iy),*)'Col. 14 - sensible heat flux,    W/m**2'
  write (n_unit(2,ix,iy),*)'Col. 15 - latent heat flux,      W/m**2'
  write (n_unit(2,ix,iy),*)'Col. 16 - downward heat flux at the upper lake surface, W/m**2'
  write (n_unit(2,ix,iy),*)'Col. 17 - momentum flux to the atmosphere, N/m**2'
  write (n_unit(2,ix,iy),*)'Col. 18 - albedo of the lake-atmosphere interface, n/d'
  
  call CHECK_UNIT(lake_series_out_unit_min,lake_series_out_unit_max, &
  & n_unit(3,ix,iy))
  open (n_unit(3,ix,iy),file=outpath(1:len_trim(outpath))//'time_series/'// &
  & 'conv_series'//coords_point//'.dat',status='unknown')
  write (n_unit(3,ix,iy),*)'Col. 1 - year'
  write (n_unit(3,ix,iy),*)'Col. 2 - month'
  write (n_unit(3,ix,iy),*)'Col. 3 - day'
  write (n_unit(3,ix,iy),*)'Col. 4 - hour'
  write (n_unit(3,ix,iy),*)'Col. 5 - the time from the start of integration, hours'
  write (n_unit(3,ix,iy),*)'Col. 6 - surface temperature, C'
  write (n_unit(3,ix,iy),*)'Col. 7 - surface density, kg/m**3'
  write (n_unit(3,ix,iy),*)'Col. 8 - heat flux downwards at the surface, K*m/s'
  write (n_unit(3,ix,iy),*)'Col. 9 - turbulent density flux at the surface, kg/m**2/s'
  write (n_unit(3,ix,iy),*)'Col. 10 - inversion depth (mass flux diagnostics),    m'
  write (n_unit(3,ix,iy),*)'Col. 11 - mixed layer depth,                          m'
  write (n_unit(3,ix,iy),*)'Col. 12 - entrainment depth,                          m'
  
! The output of LakeMIP file 1 format   
  call CHECK_UNIT(lake_series_out_unit_min,lake_series_out_unit_max, &
  & n_unit(4,ix,iy))
  open (n_unit(4,ix,iy),file=outpath(1:len_trim(outpath))//'time_series/'// &
  & 'LakeMIP_file1'//coords_point//'.dat',status='unknown')
  write (n_unit(4,ix,iy),*)'Col. 1 - time, days'
  write (n_unit(4,ix,iy),*)'Col. 2 - mixed-layer temperature &
  &(equals to the temperature of the upper water surface), degrees Celsius'
  write (n_unit(4,ix,iy),*)'Col. 3 - mean temperature of the water column, degrees Celsius'
  write (n_unit(4,ix,iy),*)'Col. 4 - bottom temperature, degrees Celsius'
  write (n_unit(4,ix,iy),*)'Col. 5 - mixed-layer depth, meters'
  write (n_unit(4,ix,iy),*)'Col. 6 - ice thickness, meters'
  write (n_unit(4,ix,iy),*)'Col. 7 - snow thickness, meters'
  write (n_unit(4,ix,iy),*)'Col. 8 - temperature at the ice upper surface, degrees Celsius'
  write (n_unit(4,ix,iy),*)'Col. 9 - temperature at the snow upper surface, degrees Celsius'
  write (n_unit(4,ix,iy),*)'Col. 10 - sensible heat flux at the lake-atmosphere &
  &interface, upwards, W/m**2'
  write (n_unit(4,ix,iy),*)'Col. 11 - latent heat flux at the lake-atmosphere &
  &interface, upwards, W/m**2'
  write (n_unit(4,ix,iy),*)'Col. 12 - momentum flux at the lake-atmosphere &
  &interface, positive, N/m**2'
  write (n_unit(4,ix,iy),*)'Col. 13 - upward long-wave radiation flux &
  &at the lake-atmosphere interface, W/m**2'
  write (n_unit(4,ix,iy),*)'Col. 14 - downward heat flux at the lake-atmosphere &
  &interface, W/m**2'
  write (n_unit(4,ix,iy),*)'Col. 15 - surface albedo, n/d'
      
  call CHECK_UNIT(lake_series_out_unit_min,lake_series_out_unit_max, &
  & n_unit(5,ix,iy))
  open (n_unit(5,ix,iy),file=outpath(1:len_trim(outpath))//'time_series/'// &
  & 'methane_series'//coords_point//'.dat',status='unknown')
  write (n_unit(5,ix,iy),*)'Col. 1 - year'
  write (n_unit(5,ix,iy),*)'Col. 2 - month'
  write (n_unit(5,ix,iy),*)'Col. 3 - day'
  write (n_unit(5,ix,iy),*)'Col. 4 - hour'
  write (n_unit(5,ix,iy),*)'Col. 5 - the time from the start of integration, hours'
  write (n_unit(5,ix,iy),*)'Col. 6 - the talik depth, m'
  write (n_unit(5,ix,iy),*)'Col. 7 - lake surface methane concentration, mol/m**3'
  write (n_unit(5,ix,iy),*)'Col. 8 - lake bottom methane concentration, mol/m**3'
  write (n_unit(5,ix,iy),*)'Col. 9 - soil bottom methane concentration, mol/m**3'
  write (n_unit(5,ix,iy),*)'Col. 10 - lake surface oxygen concentration, mol/m**3'
  write (n_unit(5,ix,iy),*)'Col. 11 - lake bottom oxygen concentration, mol/m**3'
  write (n_unit(5,ix,iy),*)'Col. 12 - total methane production due to young C decomposition, mol/(m**2*s)'
  write (n_unit(5,ix,iy),*)'Col. 13 - total methane production due to old C decomposition, mol/(m**2*s)'
  write (n_unit(5,ix,iy),*)'Col. 14 - methane ebullition flux at the lake bottom, mol/(m**2*s)'
  write (n_unit(5,ix,iy),*)'Col. 15 - methane plant-mediated flux at the lake bottom, mol/(m**2*s)'
  write (n_unit(5,ix,iy),*)'Col. 16 - methane diffusion flux at the lake bottom, upwards, mol/(m**2*s)'
  write (n_unit(5,ix,iy),*)'Col. 17 - methane turbulent flux at the lake surface, upwards, mol/(m**2*s)'
  write (n_unit(5,ix,iy),*)'Col. 18 - methane ebullition flux at the lake bottom, mg/(m**2*day)'
  write (n_unit(5,ix,iy),*)'Col. 19 - methane plant-mediated flux at the lake bottom, mg/(m**2*day)'
  write (n_unit(5,ix,iy),*)'Col. 20 - methane diffusion flux at the lake bottom, upwards, mg/(m**2*day)'
  write (n_unit(5,ix,iy),*)'Col. 21 - methane turbulent flux at the lake surface, upwards, mg/(m**2*day)'
  
  count_out(ix,iy) = int(time/(dt_out*hour_sec))
endif
      
if (int(time/(dt_out*hour_sec))>count_out(ix,iy) .or. &
  & dabs(time - (count_out(ix,iy) + 1)*dt_out*hour_sec) < small_number) then
  T_mean = varmean(Tw1(1:M+1),1)
  format_char = '(3i7,f8.2,f13.2,4f10.4,f9.2)'
  write (n_unit(1,ix,iy),format_char) year,month,day ,hour, &
  &                                   time/hour_sec ,       &
  &                                   h1, l1, hs1, ls1 
  format_char = '(3i7,f8.2,f13.2,12f11.4,f11.2)'
  write (n_unit(2,ix,iy),format_char) year,month,day,hour, &
  &                                   time/hour_sec ,      &
  &                                   tsw-273.15,Tskin(1), &   
  &                                   Tw1(1),Tw1(M+1),     &
  &                                   T_mean,T_0dim,Ti1(1),Tsn(itop), &
  &                                   hw,xlew,eflux0,tau,albedo
  format_char = '(3i7,f8.2,f13.2,2f10.4,2e15.5,3f10.4)'
  write (n_unit(3,ix,iy),format_char) year,month,day,hour,        &
  &                                   time/hour_sec,              &
  &                                   Tw1(1),row(1),eflux0_kinem, &
  &                                   turb_density_flux, zinv,    &
  &                                   H_mixed_layer, H_entrainment
  format_char = '(f10.5,10f8.2,f10.4,3f8.2)'
  write (n_unit(4,ix,iy),format_char) time/day_sec, &
  &                                   Tw1(1), T_mean, &   
  &                                   Tw1(M+1), H_mixed_layer, &
  &                                   l1, hs1, Ti1(1), &
  &                                   Tsn(itop), hw, xlew, tau, &                                 
  &                                   surfrad, eflux0, albedo
  format_char = '(3i7,f8.2,f13.2,f8.2,15e15.5)'
  write (n_unit(5,ix,iy),format_char) year,month,day,hour,        &
  &                                   time/hour_sec, h_talik,     &
  &                                   qwater(1,1), qwater(M+1,1), &
  &                                   qsoil(ns), oxyg(1,1),       &
  &                                   oxyg(M+1,1),                &
  &                                   rprod_total_newC,           &
  &                                   rprod_total_oldC,           &
  &                                   febul, fplant,              &
  &                                   - fdiff, - fdiff_lake_surf, &
  &                                   febul*mfs, fplant*mfs,      &
  &                                   - fdiff*mfs,                &
  &                                   - fdiff_lake_surf*mfs
!  format_char = '(3i7,f8.2,f13.2,f8.2,24e15.5)' ! two-meth
!  write (n_unit(5,ix,iy),format_char) year,month,day,hour,        &  ! two-meth
!  &                                   time/hour_sec, h_talik,     &  ! two-meth
!  &                                   qwater2(1,1,1:2),            &  ! two-meth
!  &                                   qwater2(M+1,1,1:2),          &  ! two-meth
!  &                                   qsoil2(ns,1:2),              &  ! two-meth
!  &                                   oxyg(1,1),                  &  ! two-meth
!  &                                   oxyg(M+1,1),                &  ! two-meth
!  &                                   rprod_total_newC,           &  ! two-meth
!  &                                   rprod_total_oldC,           &  ! two-meth
!  &                                   febul2(1:2), fplant,        &  ! two-meth
!  &                                   - fdiff2(1:2), - fdiff_lake_surf2(1:2), &  ! two-meth
!  &                                   febul2(1:2)*mfs, fplant*mfs,&  ! two-meth
!  &                                   - fdiff2(1:2)*mfs,          &  ! two-meth
!  &                                   - fdiff_lake_surf2(1:2)*mfs    ! two-meth
  
  count_out(ix,iy) = count_out(ix,iy) + 1 
endif 
      
if (firstcall(ix,iy)) firstcall(ix,iy)=.false.
END SUBROUTINE SERIES_OUT


SUBROUTINE PROFILE_OUTPUT &
& (ix, iy, nx, ny, &
&  year, month, day, hour, &
&  time, dt_out, &
&  Profile, z_prof, nprof, &
&  z_out, nout, &
&  outpath, filename, &
&  unitindex, ndec, last)

! The subroutine WATER_TEMP_OUT writes the water temperature profiles
! with the timestep dt_out

use INOUT_PARAMETERS, only : &
& lake_series_out_unit_min, &
& lake_series_out_unit_max

implicit none

! Input variables
! Reals
real(8), intent(in) :: Profile(1:nprof) ! The profile
real(8), intent(in) :: z_prof(1:nprof) ! The levels of the profile elements
real(8), intent(in) :: hour ! The hour
real(8), intent(in) :: z_out(1:nout) ! The output grid
real(8), intent(in) :: time, dt_out ! The time from beginning of integration
                                    ! The timestep of output
! Integers
integer(4), intent(in) :: ix, iy ! The current coordinates of lake point
integer(4), intent(in) :: nx, ny ! The maximal coordinates of lake points
integer(4), intent(in) :: year, month, day ! Year, month and day
integer(4), intent(in) :: nprof ! The number of layers in water, (M+1) - number of levels
integer(4), intent(in) :: nout
integer(4), intent(in) :: unitindex
integer(4), intent(in) :: ndec

character(len=*), intent(in) :: outpath
character(len=*), intent(in) :: filename

logical, intent(in) :: last

! Local variables
! Reals
real(8), parameter :: hour_sec = 60.*60.
real(8), parameter :: small_number = 1.d-5
real(8), allocatable :: Profile_out(:)


! Integers
integer(4), parameter :: unitindex_max = 20
integer(4), allocatable, save :: n_unit(:,:,:)
integer(4), allocatable, save :: count_out(:,:)
integer(4) :: i, k !Loop index

! Characters
character(len=6) :: coords_point

! Logicals
logical, allocatable, save :: firstcall(:,:)
logical :: flag


if (.not.allocated(firstcall)) then
  allocate (firstcall(1:nx, 1:ny) )
  allocate (count_out(1:nx, 1:ny) )
  allocate (n_unit   (1:unitindex_max, 1:nx, 1:ny) )
  firstcall(:,:) = .true.
  n_unit = lake_series_out_unit_min
endif  
            
if (firstcall(ix,iy)) then
  write (coords_point, '(2i3)') ix, iy

  call CHECK_UNIT(lake_series_out_unit_min,lake_series_out_unit_max, &
  & n_unit(unitindex,ix,iy))
  open (n_unit(unitindex,ix,iy), &
  & file=outpath(1:len_trim(outpath))//'time_series/'// &
  & filename(1:len_trim(filename))//coords_point//'.dat',  status='unknown')
  write (n_unit(unitindex,ix,iy),*)'Col. 1 - year'
  write (n_unit(unitindex,ix,iy),*)'Col. 2 - month'
  write (n_unit(unitindex,ix,iy),*)'Col. 3 - day'
  write (n_unit(unitindex,ix,iy),*)'Col. 4 - hour'
  write (n_unit(unitindex,ix,iy),*) &
  & 'Col. 5 - the time from the start of integration, hours'
  write (n_unit(unitindex,ix,iy),*) &
  & 'Col. [6...] - Depth1, Temp1, Depth2, Temp2, ..., DepthN, TempN'
  
  call CHECK_UNIT(lake_series_out_unit_min,lake_series_out_unit_max, &
  & n_unit(unitindex+1,ix,iy))
  open (n_unit(unitindex+1,ix,iy), &
  & file=outpath(1:len_trim(outpath))//'time_series/'// &
  & filename(1:len_trim(filename))//coords_point//'f2'//'.dat',  status='unknown')
  write (n_unit(unitindex+1,ix,iy),*)'Col. 1 - year'
  write (n_unit(unitindex+1,ix,iy),*)'Col. 2 - month'
  write (n_unit(unitindex+1,ix,iy),*)'Col. 3 - day'
  write (n_unit(unitindex+1,ix,iy),*)'Col. 4 - hour'
  write (n_unit(unitindex+1,ix,iy),*) &
  & 'Col. 5 - the time from the start of integration, hours'
  write (n_unit(unitindex+1,ix,iy),*) &
  & 'Col. 6 - Depth, m'
  write (n_unit(unitindex+1,ix,iy),*) &
  & 'Col. 7 - Temp, m'
  
  count_out(ix,iy) = int(time/(dt_out*hour_sec))
endif
      
if (int(time/(dt_out*hour_sec))>count_out(ix,iy) .or. &
  & dabs(time - (count_out(ix,iy) + 1)*dt_out*hour_sec) < small_number) then
  if (nout > 0) then
    allocate (Profile_out(1:nout))
    call LININTERPOL (z_prof,Profile,nprof,z_out,Profile_out,nout,flag)
    k = nout
    write (n_unit(unitindex,ix,iy),100) year,month,day,hour, &
    &                                   time/hour_sec,       &
    &                                   ((z_out(i),Profile_out(i)), i = 1, nout)
    do i = 1, nout
      write (n_unit(unitindex+1,ix,iy),101) year,month,day,hour, &
      &                                   time/hour_sec,         &
      &                                   -z_out(i),Profile_out(i)
    enddo
    deallocate (Profile_out)
  else
    k = nprof
    write (n_unit(unitindex,ix,iy),100) year,month,day,hour, &
    &                                   time/hour_sec,       &
    &                                   ((z_prof(i),Profile(i)), i = 1, nprof)
    do i = 1, nprof
      write (n_unit(unitindex+1,ix,iy),101) year,month,day,hour, &
      &                                   time/hour_sec,         &
      &                                   -z_prof(i),Profile(i)
    enddo
  endif
  if (last) count_out(ix,iy) = count_out(ix,iy) + 1 
endif 
      
if (firstcall(ix,iy) .and. last) firstcall(ix,iy) = .false.
100 format (3i7,f8.2,f13.2,<k>(f10.2, f<6+ndec>.<ndec>))
101 format (3i7,f8.2,f13.2,f10.2,f<6+ndec>.<ndec>)
END SUBROUTINE PROFILE_OUTPUT


SUBROUTINE ACCUM_VAR &
& (dt, l1, febul, fdiff, fdiff_lake_surf, &
& rprod_total_newC, rprod_total_oldC, &
& febultot, fdifftot, fdiff_lake_surftot, &
& rprod_total_newC_integr, rprod_total_oldC_integr, &
& add_to_winter)

use METH_OXYG_CONSTANTS, only : &
& molmass_ch4

implicit none

! Input/output variables

real(8), intent(in) :: dt
real(8), intent(in) :: l1
real(8), intent(in) :: febul
real(8), intent(in) :: fdiff
real(8), intent(in) :: fdiff_lake_surf
!real(8), intent(in) :: febul(1:2) ! two-meth
!real(8), intent(in) :: fdiff(1:2) ! two-meth
!real(8), intent(in) :: fdiff_lake_surf(1:2) ! two-meth
real(8), intent(in) :: rprod_total_newC
real(8), intent(in) :: rprod_total_oldC
logical, intent(inout) :: add_to_winter
real(8), intent(inout) :: febultot(1:2)
real(8), intent(inout) :: fdifftot(1:2)
real(8), intent(inout) :: fdiff_lake_surftot(1:2)
real(8), intent(inout) :: rprod_total_newC_integr(1:2)
real(8), intent(inout) :: rprod_total_oldC_integr(1:2)
!real(8), intent(inout) :: febultot(1:2,1:2) ! two-meth
!real(8), intent(inout) :: fdifftot(1:2,1:2) ! two-meth
!real(8), intent(inout) :: fdiff_lake_surftot(1:2,1:2) ! two-meth
!real(8), intent(inout) :: rprod_total_newC_integr(1:2,1:2) ! two-meth
!real(8), intent(inout) :: rprod_total_oldC_integr(1:2,1:2) ! two-meth

! Local variables

real(8), parameter :: mfs = molmass_ch4*8.64d+7 ! the transform multiplier for methane flux
                                                ! from mol/(m**2*s) to mg/(m**2*day)
real(8), parameter :: day_sec = 24.*60.*60.

integer(4) :: i ! loop index

if (l1 == 0. .and. (.not. add_to_winter)) then
  febultot(1) = febultot(1) + febul*mfs*dt/day_sec
  fdifftot(1) = fdifftot(1) - fdiff*mfs*dt/day_sec
  fdiff_lake_surftot(1) = fdiff_lake_surftot(1) - fdiff_lake_surf*mfs*dt/day_sec
  rprod_total_newC_integr(1) = rprod_total_newC_integr(1) + rprod_total_newC*mfs*dt/day_sec 
  rprod_total_oldC_integr(1) = rprod_total_oldC_integr(1) + rprod_total_oldC*mfs*dt/day_sec 
else
  febultot(2) = febultot(2) + febul*mfs*dt/day_sec
  fdifftot(2) = fdifftot(2) - fdiff*mfs*dt/day_sec
  fdiff_lake_surftot(2) = fdiff_lake_surftot(2) - fdiff_lake_surf*mfs*dt/day_sec
  rprod_total_newC_integr(2) = rprod_total_newC_integr(2) + rprod_total_newC*mfs*dt/day_sec 
  rprod_total_oldC_integr(2) = rprod_total_oldC_integr(2) + rprod_total_oldC*mfs*dt/day_sec
  if (add_to_winter) add_to_winter = .false.
endif

!do i = 1, 2 ! two-meth
!  if (l1 == 0.) then ! two-meth
!    febultot(1,i) = febultot(1,i) + febul(i)*mfs*dt/day_sec ! two-meth
!    fdifftot(1,i) = fdifftot(1,i) - fdiff(i)*mfs*dt/day_sec ! two-meth
!    fdiff_lake_surftot(1,i) = fdiff_lake_surftot(1,i) - fdiff_lake_surf(i)*mfs*dt/day_sec ! two-meth
!    rprod_total_newC_integr(1,i) = rprod_total_newC_integr(1,i) + rprod_total_newC*mfs*dt/day_sec  ! two-meth
!    rprod_total_oldC_integr(1,i) = rprod_total_oldC_integr(1,i) + rprod_total_oldC*mfs*dt/day_sec  ! two-meth
!  else ! two-meth
!    febultot(2,i) = febultot(2,i) + febul(i)*mfs*dt/day_sec ! two-meth
!    fdifftot(2,i) = fdifftot(2,i) - fdiff(i)*mfs*dt/day_sec ! two-meth
!    fdiff_lake_surftot(2,i) = fdiff_lake_surftot(2,i) - fdiff_lake_surf(i)*mfs*dt/day_sec ! two-meth
!    rprod_total_newC_integr(2,i) = rprod_total_newC_integr(2,i) + rprod_total_newC*mfs*dt/day_sec  ! two-meth
!    rprod_total_oldC_integr(2,i) = rprod_total_oldC_integr(2,i) + rprod_total_oldC*mfs*dt/day_sec  ! two-meth
!  endif ! two-meth
!enddo ! two-meth

END SUBROUTINE ACCUM_VAR    


SUBROUTINE CHECK_UNIT(unit_min,unit_max,nunit)
! The subroutine CHECK_UNIT checks if the output unit is already occupied,
! and if yes, it returns the number of the free unit
implicit none

! Input variables
integer(4), intent(in) :: unit_min
integer(4), intent(in) :: unit_max

! Input/output variables
integer(4), intent(inout) :: nunit

! Local variables
logical :: unit_opened

inquire (unit=nunit,opened=unit_opened)

if (unit_opened) then
  do while (unit_opened)
!    write(*,*) 'The unit ', nunit, 'is attempted &
!    & to be connected to a file, while already connected: incrementing unit'
    nunit = nunit + 1
    inquire (unit=nunit,opened=unit_opened)
  enddo
!  STOP
endif

if (nunit < unit_min .or. nunit > unit_max) then
  write(*,*) 'Error on LAKE model: the bounds of permitted input/ouput &
  &unit numbers exceeded: STOP'
  STOP
endif

END SUBROUTINE CHECK_UNIT


FUNCTION GETVARVAL(n1,n2,line,name)
implicit none

real(8) :: GETVARVAL

! Input variables
integer,       intent(in):: n2,n1
character*200, intent(in):: line
character*25,  intent(in):: name

! Local variables
real(8) work

read (line((n2+1):100),*) work
!print*, name//' = ', work

GETVARVAL = work      

RETURN
END FUNCTION GETVARVAL


FUNCTION IGETVARVAL(n1,n2,line,name)
implicit none

integer(4) :: IGETVARVAL

! Input variables
integer(4),    intent(in):: n2,n1
character*200, intent(in):: line
character*25,  intent(in):: name

! Local variables
integer(4) iwork

read (line((n2+1):100),*) iwork
!print*, name//' = ', iwork

IGETVARVAL = iwork      

RETURN
END FUNCTION IGETVARVAL
